{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd5d45c3-c7f5-4eb0-b1a2-74e91e4a54e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import struct\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.path import Path\n",
    "import matplotlib as mpl\n",
    "import re\n",
    "#import fiona\n",
    "import cartopy.crs as ccrs\n",
    "import matplotlib.ticker as mticker\n",
    "from scipy import interpolate\n",
    "\n",
    "\n",
    "%config InlineBackend.figure_format='retina'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "643be159-5631-4749-81fc-f5d661b0bd9c",
   "metadata": {},
   "source": [
    "## Load Datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e44b570-7002-4bcf-a4a4-560e405d99c8",
   "metadata": {
    "tags": []
   },
   "source": [
    "### LROC WAC Mosaic\n",
    "\n",
    "Note: using PDS3 version but recently has been updated to PDS4. Have not yet figure out how to display .TIF file.\n",
    "\n",
    "https://pds.lroc.asu.edu/data/LRO-L-LROC-5-RDR-V1.0/LROLRC_2001/EXTRAS/BROWSE/WAC_GLOBAL/WAC_GLOBAL_E000N0000_004P.TIF\n",
    "\n",
    "First, remove header lines and save to new file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "713948de-22fd-4ea3-88ac-7c6a970ca28b",
   "metadata": {},
   "outputs": [],
   "source": [
    "wacfile = '/work1/moon/WAC/WAC_GLOBAL_E000N0000_004P_nohead.txt'\n",
    "\n",
    "\n",
    "LINES = 720 #rows / lat\n",
    "LINE_SAMPLES = 1440 #collumns /lon\n",
    "SAMPLE_BITS = 32\n",
    "BYTES_PER_PIXEL = int(SAMPLE_BITS/8)\n",
    "\n",
    "llat = np.linspace(-90,90,LINES)\n",
    "llon = np.linspace(-180,180,LINE_SAMPLES)\n",
    "LLON, LLAT = np.meshgrid(llon,llat)\n",
    "\n",
    "pixels = LINES*LINE_SAMPLES\n",
    "fmt = '{endian}{pixels}{fmt}'.format(endian='<',pixels=pixels,fmt='f')\n",
    "#plt.figure(figsize = [20,16])\n",
    "with open(wacfile,'rb') as f:\n",
    "    wacimg = np.array(struct.unpack(fmt,f.read(pixels*BYTES_PER_PIXEL))).reshape(LINES,LINE_SAMPLES)\n",
    "wacimg_corr = np.flip(np.flip(wacimg), axis=1)\n",
    "wacimg_corr_adj = np.roll(wacimg_corr, 1440-210,axis =1) ## manual adjustment for strange offset - this is for display only!!!!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "586c465e-dc48-4083-bf41-0f8539c6f780",
   "metadata": {},
   "source": [
    "### Imbrium Ejecta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c390d12b-b811-4695-a77b-1d68c5422a76",
   "metadata": {},
   "outputs": [],
   "source": [
    "ej_th = np.genfromtxt('../results/imbrium721_45.txt',skip_header=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "221da0fa-161e-47b3-a395-de4c83545709",
   "metadata": {
    "tags": []
   },
   "source": [
    "### TSTA Results File:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b0ba328-3db1-4c1b-9714-b9d445c4b12a",
   "metadata": {},
   "outputs": [],
   "source": [
    "full_craterdat = pd.read_csv('../results/tsta_20_200.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdc8f232-f4f4-4783-9052-49e59e1f623d",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Constants and Functions\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d8fb2a5-1df5-437d-a7ca-d9f98655f9e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "## function to sort points counter clockwise (for plotting)\n",
    "def ccw_sort(p):\n",
    "    p = np.array(p)\n",
    "    mean = np.mean(p,axis=0)\n",
    "    d = p-mean\n",
    "    s = np.arctan2(d[:,0], d[:,1])\n",
    "    return p[np.argsort(s),:]\n",
    "\n",
    "## function for great-circle distance between points\n",
    "def great_circle_dist_moon(p1,p2,r=1737.0):\n",
    "    \"\"\"\n",
    "    Calculate great-circle distance between two points\n",
    "    \n",
    "    Input:\n",
    "        p1 = [lon, lat] \n",
    "        p2 = [lon, lat]\n",
    "\n",
    "        r = Radius (Default 1737.0 km) \n",
    "    Output:\n",
    "        d = distance (Default km)\n",
    "    \"\"\"\n",
    "    \n",
    "    lon1 = np.deg2rad(p1[0])\n",
    "    lon2 = np.deg2rad(p2[0])\n",
    "    lat1 = np.deg2rad(p1[1])\n",
    "    lat2 = np.deg2rad(p2[1])\n",
    "    \n",
    "    cent_ang = np.arccos(np.sin(lat1)*np.sin(lat2) + np.cos(lat1)*np.cos(lat2)*np.cos(np.abs(lon1-lon2)))\n",
    "    d = r*cent_ang\n",
    "    return d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab5bbdae-4fcc-4c8f-a085-cdead4619a3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "th_c = 0.63 #ppm metzger 1973\n",
    "\n",
    "rho_c = 2550.0 #kg/m**3 Wieczorek 2013  bulk density of highlands crust\n",
    "rho_c_km = rho_c*(1000**3) #kg/km**3\n",
    "\n",
    "D_sc = 18 #km (simple complex transition)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "671de2d9-758b-46ca-bffb-8b5687f34ae3",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Imbrium (from Neumann 2015)\n",
    "imb_lat = 37.0\n",
    "imb_lon = 341.5 - 360\n",
    "imb_dc = 1321 # main ring \n",
    "dist_from_imb = great_circle_dist_moon([imb_lon,imb_lat],[LLON,LLAT])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cab7ca74-4cfa-4cd4-85e5-d389f75fc895",
   "metadata": {},
   "outputs": [],
   "source": [
    "####functions for spherical coords\n",
    "def vec_rep(lon,lat):\n",
    "    \"\"\"\n",
    "    Cartesian vector on unit sphere from lat,lon\n",
    "    Multiply output by radius of body if necessary\n",
    "    \"\"\"\n",
    "    phi = np.deg2rad(lon)\n",
    "    thet = np.deg2rad(90. -lat) ##colatitude\n",
    "    \n",
    "    #point = np.zeros((len(lon),3)) ##cartesian vector\n",
    "    x = np.sin(thet)*np.cos(phi)\n",
    "    y = np.sin(thet)*np.sin(phi)\n",
    "    z = np.cos(thet)\n",
    "    return x,y,z\n",
    "\n",
    "def rodr_rot(v,theta,k):\n",
    "    \"\"\"\n",
    "    Rotate vector v around unit vector k (representing axis of rotation) by theta degrees\n",
    "  \n",
    "    \"\"\"\n",
    "    \n",
    "    v_rot = v*np.cos(theta) + np.sin(theta)*np.cross(k,v) + np.outer(np.dot(v,k)*(1-np.cos(theta)),k)\n",
    "    return v_rot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2220a566-07a2-4a60-add6-4b9b5bb183b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "###scenario 1\n",
    "def kreep_over(dat,quad,th_k,X):\n",
    "    \"\"\"\n",
    "    dat = pandas dataset \n",
    "    quad = subset of pandas dataset\n",
    "    th_k = numpy array KREEP thorium concentration\n",
    "    X = numpy array Thickness of KREEP layer\n",
    "    \"\"\"\n",
    "    XX,tt = np.meshgrid(X,th_k)\n",
    "\n",
    "    data = np.zeros((len(th_k),len(X),len(dat.loc[quad])))\n",
    "    error = np.zeros_like(data)\n",
    "    a = 0\n",
    "    for _,rows in dat.loc[quad].iterrows():\n",
    "        x_temp = np.zeros_like(XX) #array height of small paraboloid \n",
    "        vc_temp = np.zeros_like(XX)\n",
    "        vk_temp = np.zeros_like(XX)\n",
    "        x_temp[np.where(XX < rows['h_exc'])] = rows['h_exc'] - XX[np.where(XX < rows['h_exc'])] #array height of small paraboloid \n",
    "        vc_temp = 5*np.pi*rows['r_at']*x_temp**2/2.6 ##volume of small paraboloid\n",
    "        vk_temp = -vc_temp + rows['vols']\n",
    "        data[:,:,a] = (tt*vk_temp + th_c*vc_temp)/rows['vols']\n",
    "        error[:,:,a] = data[:,:,a] - rows['ejecta_m_t']\n",
    "        a += 1\n",
    "\n",
    "\n",
    "    rms_err = np.sqrt(np.mean(error**2,axis = 2))\n",
    "    \n",
    "    return rms_err\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "481316d2-37d4-4788-83c5-7b056491752b",
   "metadata": {},
   "outputs": [],
   "source": [
    "###vertical axis (imb-centric) \n",
    "\n",
    "vert_lat = np.arange(-90.,90.1,.5)\n",
    "vert_lon = np.zeros_like(vert_lat) + imb_lon\n",
    "\n",
    "#horizontal axis\n",
    "lon_sp = np.arange(-180.,180.1,1)\n",
    "ha_cx,ha_cy,ha_cz = vec_rep(lon_sp,np.zeros_like(lon_sp))\n",
    "\n",
    "ha_d = np.zeros((len(ha_cx),3))\n",
    "ha_d[:,0] = ha_cx\n",
    "ha_d[:,1] = ha_cy\n",
    "ha_d[:,2] = ha_cz\n",
    "\n",
    "\n",
    "##undo rotation \n",
    "\n",
    "j = np.array([0.,1.,0.])\n",
    "undo_j = rodr_rot(ha_d,-np.deg2rad(np.abs(imb_lat)),j)\n",
    "\n",
    "k = np.array([0.,0.,1.])\n",
    "undo_k = rodr_rot(undo_j, -np.deg2rad(np.abs(imb_lon)),k)\n",
    "\n",
    "ha_thet = np.arccos(undo_k[:,2]/np.sqrt(undo_k[:,0]**2 + undo_k[:,1]**2 + undo_k[:,2]**2))\n",
    "\n",
    "ha_phi = np.sign(undo_k[:,1])*np.arccos(undo_k[:,0]/np.sqrt(undo_k[:,0]**2 + undo_k[:,1]**2))\n",
    "\n",
    "ha_lat = 90. - np.rad2deg(ha_thet)\n",
    "ha_lon = np.rad2deg(ha_phi)\n",
    "\n",
    "###latitude axis\n",
    "h_idx = np.argsort(ha_lon)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3e2c976-e41f-49bf-9950-0aa6a800cb0f",
   "metadata": {},
   "source": [
    "### Imbrium-centric analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f4ae01c-8366-406f-8188-776bc48d7412",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "full_craterdat = full_craterdat.loc[full_craterdat['floor_m_t'].notna()& full_craterdat['ejecta_m_t'].notna()] #exclude nan locations\n",
    "\n",
    "##distance from imbrium\n",
    "full_craterdat['dist_from_imbrium'] = great_circle_dist_moon([imb_lon,imb_lat],[full_craterdat['lon'],full_craterdat['lat']])\n",
    "\n",
    "##apparent transient crater\n",
    "full_craterdat['r_at'] = (full_craterdat['rc']**(0.921)*(D_sc/2)**(0.079))/(1.32) #adapted from holsapple 1993 with Rat = Rt/1.3\n",
    "\n",
    "#excavation depth\n",
    "full_craterdat['h_exc'] = 1.3*full_craterdat['r_at']/5.0 #melosh pg. 80\n",
    "\n",
    "full_craterdat['diff'] = full_craterdat['ejecta_m_t'] - full_craterdat['floor_m_t'] \n",
    "full_craterdat['vols'] = np.pi*full_craterdat['h_exc']*full_craterdat['r_at']**2/2\n",
    "\n",
    "\n",
    "\n",
    "full_craterdat = full_craterdat.loc[(full_craterdat['ejecta_m_t'] > 0) & (full_craterdat['floor_m_t'] > 0)]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd94025e-598f-4418-8489-0154daefad10",
   "metadata": {},
   "outputs": [],
   "source": [
    "###rotating the reference frame \n",
    "\n",
    "### convert to cartesian\n",
    "\n",
    "c_x,c_y,c_z = vec_rep(full_craterdat[\"lon\"].values, full_craterdat[\"lat\"].values)\n",
    "\n",
    "cart_data = np.zeros((len(c_x),3))\n",
    "cart_data[:,0] = c_x\n",
    "cart_data[:,1] = c_y\n",
    "cart_data[:,2] = c_z\n",
    "\n",
    "### first need to rotate around z by lon\n",
    "k = np.array([0.,0.,1.])\n",
    "rotz = rodr_rot(cart_data,np.deg2rad(np.abs(imb_lon)),k)\n",
    "\n",
    "### now rotate around y by lat\n",
    "j = np.array([0.,1.,0.])\n",
    "rotzy = rodr_rot(rotz,np.deg2rad(np.abs(imb_lat)),j)\n",
    "i_thet = np.arccos(rotzy[:,2]/np.sqrt(rotzy[:,0]**2 + rotzy[:,1]**2 + rotzy[:,2]**2))\n",
    "\n",
    "i_phi = np.sign(rotzy[:,1])*np.arccos(rotzy[:,0]/np.sqrt(rotzy[:,0]**2 + rotzy[:,1]**2))\n",
    "i_lat = 90. - np.rad2deg(i_thet)\n",
    "i_lon = np.rad2deg(i_phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b02278c6-bf99-40ba-b72a-83ff77b5cad9",
   "metadata": {},
   "outputs": [],
   "source": [
    "imb_quad = np.zeros_like(i_lat)\n",
    "\n",
    "imb_quad[np.where((i_lat > 0.) & (i_lon > 0.))] = 1\n",
    "imb_quad[np.where((i_lat > 0.) & (i_lon < 0.))] = 2\n",
    "imb_quad[np.where((i_lat < 0.) & (i_lon < 0.))] = 3\n",
    "imb_quad[np.where((i_lat < 0.) & (i_lon > 0.))] = 4\n",
    "\n",
    "\n",
    "#full_craterdat['i_lon'] = ic_lon\n",
    "#full_craterdat['i_lat'] = ic_lat\n",
    "full_craterdat['imb_quad'] = imb_quad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "968172d7-2a32-4d1d-849a-65583be3c748",
   "metadata": {},
   "outputs": [],
   "source": [
    "contours = np.arange(0,5000,400)\n",
    "plt.pcolormesh(LLON,LLAT,wacimg_corr_adj, cmap = 'gray', vmin = -0.03, vmax = 0.2, shading='auto')\n",
    "plt.plot(vert_lon,vert_lat)\n",
    "plt.plot(ha_lon[h_idx],ha_lat[h_idx])\n",
    "\n",
    "colors = np.array(['red','orange','yellow','green'])\n",
    "\n",
    "color_idx = np.array(full_craterdat['imb_quad'].values - 1,dtype = int)\n",
    "plt.scatter(full_craterdat['lon'],full_craterdat['lat'], color = colors[color_idx])\n",
    "\n",
    "c = plt.contour(LLON,LLAT,dist_from_imb,contours,colors = 'white' )\n",
    "plt.clabel(c, inline=1, fontsize=12,fmt='%d')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ed02cf7-89b4-4707-a2f4-25e382eaf19c",
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_list = np.arange(0,5000,500) \n",
    "\n",
    "fig = plt.figure(figsize = [9,4])\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "ax.set_yscale('log')\n",
    "ax.hist([full_craterdat['dist_from_imbrium'].loc[full_craterdat['diff'] > 0],\n",
    "          full_craterdat['dist_from_imbrium'].loc[full_craterdat['diff'] < 0]],bins = bin_list,\n",
    "         color = ['red','blue'],alpha = 0.2,label = ['Total S1', 'Total S2'])\n",
    "ax.hist([full_craterdat['dist_from_imbrium'].loc[full_craterdat['diff'] > full_craterdat['back_std']],\n",
    "          full_craterdat['dist_from_imbrium'].loc[full_craterdat['diff'] < -full_craterdat['back_std']]],bins = bin_list,\n",
    "         color = ['red','blue'],label = ['Significant S1', 'Significant S2'])\n",
    "\n",
    "ax.plot([imb_dc/2,imb_dc/2],[0,70],'k--',linewidth = 1,label = \"Imbrium Main Ring\")\n",
    "ax.set_xlabel('Distance from center of Imbrium (km)')\n",
    "ax.set_ylabel('Number of Craters')\n",
    "ax.legend(loc = 'upper left')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17fcf7bf-1f5a-45c5-9e7f-ea2eb08d5d9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.savefig('craters_w_dist_from_imb_tot.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9f4a3c0-6a8f-45c1-b3a7-46d4e04dfee0",
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_list = np.arange(0,5000,500) \n",
    "fig,ax = plt.subplots(nrows = 4,figsize = [9,8])\n",
    "for i in range(4):\n",
    "    ax[i].set_yscale('log')\n",
    "    quad = full_craterdat.loc[full_craterdat['imb_quad'] == i + 1]\n",
    "    ax[i].hist([quad['dist_from_imbrium'].loc[quad['diff'] > 0],\n",
    "          quad['dist_from_imbrium'].loc[quad['diff'] < 0]],bins = bin_list,\n",
    "         color = ['red','blue'],alpha = 0.2,label = ['Total S1', 'Total S2'])\n",
    "    ax[i].hist([quad['dist_from_imbrium'].loc[quad['diff'] > quad['back_std']],\n",
    "          quad['dist_from_imbrium'].loc[quad['diff'] < -quad['back_std']]],bins = bin_list,\n",
    "         color = ['red','blue'],label = ['Significant S1', 'Significant S2'])\n",
    "    if i < 3:\n",
    "        ax[i].set_xticks([])\n",
    "    else:\n",
    "        ax[i].set_xlabel('Distance from center of Imbrium (km)')\n",
    "    ax[i].set_ylabel('Craters in Q{}'.format(i+1))\n",
    "#ax.plot([imb_dc/2,imb_dc/2],[0,70],'k--',linewidth = 1,label = \"Imbrium Main Ring\")\n",
    "#ax.set_xlabel('Distance from center of Imbrium (km)')\n",
    "#ax.set_ylabel('Number of Craters')\n",
    "#ax.legend(loc = 'upper left')\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f74973ac-f880-439a-b41c-ae88f3f6e644",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Imbrium mass balance\n",
    "th_k = np.linspace(th_c,30,150) # KREEP thorium concentration\n",
    "\n",
    "dist_list = np.arange(500,5000,200)\n",
    "\n",
    "depths = np.zeros((4,len(dist_list)))\n",
    "thors = np.zeros((4,len(dist_list)))\n",
    "for q in range(4):\n",
    "    quad = full_craterdat.loc[full_craterdat['imb_quad'] == q + 1]\n",
    "    \n",
    "    i = 0\n",
    "    for im_lim in dist_list:\n",
    "        imb_b = [im_lim - 100,im_lim + 100]\n",
    "        reg_s1 = (quad['dist_from_imbrium'] > imb_b[0]) & (quad['dist_from_imbrium'] < imb_b[1]) &\\\n",
    "        (quad['diff'] >= quad['back_std'])\n",
    "\n",
    "        max_h = quad.loc[reg_s1]['h_exc'].max()\n",
    "    \n",
    "    \n",
    "        X = np.linspace(0,max_h,1000) # Thickness of KREEP layer\n",
    "        XX,tt = np.meshgrid(X,th_k)\n",
    "\n",
    "\n",
    "        rms_err_r = kreep_over(quad,reg_s1,th_k,X)\n",
    "\n",
    "        max_ejr = quad.loc[reg_s1]['ejecta_m_t'].max()\n",
    "\n",
    "        max_ej_idxr = np.argmin(np.abs(th_k - max_ejr))\n",
    "\n",
    "        max_ej_errors = rms_err_r[max_ej_idxr,:]\n",
    "        max_ej_idx = np.argmin(max_ej_errors)\n",
    "        max_thick = X[max_ej_idx]\n",
    "\n",
    "    \n",
    "        depths[q,i] = max_thick\n",
    "        if np.isnan(max_thick):\n",
    "            thors[q,i] = np.nan\n",
    "        else:\n",
    "            thors[q,i] = max_ejr\n",
    "        i +=1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "178708ce-9cdc-49f3-b4aa-b0856fd99510",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(nrows = 2, ncols =2, figsize = [18,12])\n",
    "\n",
    "for i in range(4):\n",
    "    ax[i//2,i%2].plot(dist_list,depths[i,:],'k*')\n",
    "    ax[i//2,i%2].set_title('Quadrant {}'.format(i + 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23135c83-f177-4b39-97cf-5fbb563bbec0",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.zeros((3,3)) + np.array([1,2,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17f70eca-298d-4360-aad1-54962576215e",
   "metadata": {},
   "outputs": [],
   "source": [
    "####SSE of fit for a bunch of different model parameters \n",
    "kreep_val = np.arange(3,70,.5)\n",
    "\n",
    "bounds = [imb_dc,3300]\n",
    "an_th = 0.63  ###from metzger 1973\n",
    "\n",
    "sse_tot = np.zeros((1,len(kreep_val)))\n",
    "sse_qs = np.zeros((4,len(kreep_val)))\n",
    "\n",
    "idcs = np.where((dist_list < bounds[1])& (dist_list > bounds[0]))\n",
    "\n",
    "for i in range(len(kreep_val)):\n",
    "    kr_th = kreep_val[i]\n",
    "    \n",
    "    con = (ej_th[:,2]*kr_th + ej_th[:,3]*an_th)/ej_th[:,1] \n",
    "    s = interpolate.InterpolatedUnivariateSpline(ej_th[:,0], con) ###interpolate the data points with a univariate spline\n",
    "    \n",
    "    \n",
    "    er_qs = np.zeros((4,len(idcs))) + s(dist_list[idcs])\n",
    "    \n",
    "    for q in range(4):\n",
    "        er_qs[q,:] += -thors[q,idcs][0] \n",
    "    \n",
    "    sse_qs[:,i] = np.nansum(er_qs**2,axis = 1)\n",
    "\n",
    "    sse_tot[0,i] = np.nansum(er_qs**2)\n",
    "ej_inc = kreep_val[np.argmin(sse_tot)]\n",
    "\n",
    "####minimum error w/ bounds \n",
    "\n",
    "min_err = np.min(sse_tot[0,:])\n",
    "\n",
    "###+ 100% error bounds\n",
    "err100 = np.sort(kreep_val[np.argsort(np.abs(sse_tot[0,:] - 2*min_err))[:2]])\n",
    "err200 = np.sort(kreep_val[np.argsort(np.abs(sse_tot[0,:] - 4*min_err))[:2]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f52c6db-6e70-44eb-a788-ab77997ec49d",
   "metadata": {},
   "outputs": [],
   "source": [
    "err200"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52791df8-a5d8-4d9f-a50f-423afe10bf0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "commb_no = (ej_th[:,2]*ej_inc + ej_th[:,3]*0.63)/ej_th[:,1] \n",
    "err200_1 = (ej_th[:,2]*err200[0] + ej_th[:,3]*0.63)/ej_th[:,1] \n",
    "err200_2 = (ej_th[:,2]*err200[1] + ej_th[:,3]*0.63)/ej_th[:,1] \n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(nrows = 2, figsize = [8,10])\n",
    "ax[0].semilogy([3300,3300],[0,15], color = 'blue', linestyle = '--',alpha = 0.6)\n",
    "ax[0].text(3320,0.025,'Minimum Modeled Thickness',rotation=90,color = 'blue', alpha = 0.6)\n",
    "ax[0].semilogy(ej_th[:,0],ej_th[:,1]/1000,'k-', linewidth = 3, label = 'Total Deposits' )\n",
    "ax[0].semilogy(ej_th[:,0],ej_th[:,3]/1000, linestyle = '--',color = 'dimgray',linewidth = 3, label = 'Local Material')\n",
    "ax[0].semilogy(ej_th[:,0],ej_th[:,2]/1000,linestyle = ':', color = 'darkgrey',linewidth = 3, label = 'Primary Ejecta')\n",
    "\n",
    "colors = ['maroon','orangered','darkorange','goldenrod']\n",
    "markers = ['*','o','<','>']\n",
    "for q in range(4):\n",
    "    ax[0].semilogy(dist_list,depths[q,:],marker = markers[q], color = colors[q],linestyle = 'None', label = 'Quadrant {}'.format(q + 1), markersize = 10)\n",
    "    ax[1].plot(dist_list,thors[q,:],marker = markers[q], color = colors[q],linestyle = 'None',label = 'Quadrant {}'.format(q + 1),markersize = 10)\n",
    "\n",
    "ax[0].set_ylim([0.02,15])\n",
    "ax[0].set_xlim([500,5100])\n",
    "\n",
    "ax[0].set_ylabel('Thickness (km)', fontsize = 12)\n",
    "ax[0].legend(loc='lower left', fontsize = 14)\n",
    "ax[0].tick_params(axis='both', which='major', labelsize=12)\n",
    "\n",
    "\n",
    "ax[1].plot([imb_dc,imb_dc],[0,25], color = 'blue', linestyle = '--',alpha = 0.6)\n",
    "ax[1].plot([3300,3300],[0,25], color = 'blue', linestyle = '--',alpha = 0.6)\n",
    "#ax[1].text(3330,4,'Minimum  Thickness',rotation=90,color = 'blue', alpha = 0.4)\n",
    "#ax[1].text(imb_dc + 10,10.,'Approx. End of Continuous Ejecta',rotation=90,color = 'blue', alpha = 0.6)\n",
    "\n",
    "ax[1].plot(ej_th[:,0],commb_no,'k-',linewidth= 3,label = 'Best Fit: {:.0f} ppm'.format(ej_inc))\n",
    "ax[1].fill_between(ej_th[:,0],err200_1,err200_2,color = 'grey',alpha = 0.2,\n",
    "             label = '$\\pm$200% error: \\n {:.0f} ppm - {:.0f} ppm'.format(err200[0], err200[1]))\n",
    "\n",
    "\n",
    "ax[1].set_xlabel('Distance from Center of Imbrium', fontsize = 14)\n",
    "ax[1].set_xlim([500,5100])\n",
    "ax[1].legend(loc='upper right', fontsize = 14)\n",
    "\n",
    "\n",
    "ax[1].set_ylabel('Thorium (ppm)',fontsize = 14)\n",
    "ax[1].tick_params(axis='both', which='major', labelsize=12)\n",
    "\n",
    "\n",
    "\n",
    "ax[1].set_ylim([0,22])\n",
    "#ax[1].legend(loc='upper right',fontsize = 14)\n",
    "#ax[1].set_title('Minimum Concentration of Layer',fontsize = 14)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bdc270a-4b22-4ad7-921d-8fea41744a6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.savefig('imbrium_figure_corrected_2024.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc76546b-8b47-42bf-b131-011d457a6d8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "commb_25 = (ej_th[:,2]*25 + ej_th[:,3]*0.63)/ej_th[:,1]\n",
    "commb_35 = (ej_th[:,2]*35 + ej_th[:,3]*0.63)/ej_th[:,1]\n",
    "commb_13 = (ej_th[:,2]*13 + ej_th[:,3]*0.63)/ej_th[:,1]\n",
    "\n",
    " \n",
    "err200_1 = (ej_th[:,2]*err200[0] + ej_th[:,3]*0.63)/ej_th[:,1] \n",
    "err200_2 = (ej_th[:,2]*err200[1] + ej_th[:,3]*0.63)/ej_th[:,1] \n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(nrows = 2, figsize = [8,10])\n",
    "ax[0].semilogy([3300,3300],[0,15], color = 'blue', linestyle = '--',alpha = 0.6)\n",
    "ax[0].text(3320,0.025,'Minimum Modeled Thickness',rotation=90,color = 'blue', alpha = 0.6)\n",
    "ax[0].semilogy(ej_th[:,0],ej_th[:,1]/1000,'k-', linewidth = 3, label = 'Total Deposits' )\n",
    "ax[0].semilogy(ej_th[:,0],ej_th[:,3]/1000, linestyle = '--',color = 'dimgray',linewidth = 3, label = 'Local Material')\n",
    "ax[0].semilogy(ej_th[:,0],ej_th[:,2]/1000,linestyle = ':', color = 'darkgrey',linewidth = 3, label = 'Primary Ejecta')\n",
    "\n",
    "colors = ['maroon','orangered','darkorange','goldenrod']\n",
    "markers = ['*','o','<','>']\n",
    "for q in range(4):\n",
    "    ax[0].semilogy(dist_list,depths[q,:],marker = markers[q], color = colors[q],linestyle = 'None', label = 'Quadrant {}'.format(q + 1), markersize = 10)\n",
    "    ax[1].plot(dist_list,thors[q,:],marker = markers[q], color = colors[q],linestyle = 'None',label = 'Quadrant {}'.format(q + 1),markersize = 10)\n",
    "\n",
    "ax[0].set_ylim([0.02,15])\n",
    "ax[0].set_xlim([500,5100])\n",
    "\n",
    "ax[0].set_ylabel('Thickness (km)', fontsize = 12)\n",
    "ax[0].legend(loc='lower left', fontsize = 14)\n",
    "ax[0].tick_params(axis='both', which='major', labelsize=12)\n",
    "\n",
    "\n",
    "ax[1].plot([imb_dc,imb_dc],[0,25], color = 'blue', linestyle = '--',alpha = 0.6)\n",
    "ax[1].plot([3300,3300],[0,25], color = 'blue', linestyle = '--',alpha = 0.6)\n",
    "#ax[1].text(3330,4,'Minimum  Thickness',rotation=90,color = 'blue', alpha = 0.4)\n",
    "ax[1].text(imb_dc + 10,10.,'Edge of Continuous Ejecta',rotation=90,color = 'blue', alpha = 0.6)\n",
    "\n",
    "#ax[1].plot(ej_th[:,0],commb_25,'k-',linewidth= 3,label = 'Best Fit: {:.0f} ppm'.format(ej_inc))\n",
    "#ax[1].plot(ej_th[:,0],commb_35,'k-',linewidth= 3,label = 'Best Fit: {:.0f} ppm'.format(ej_inc))\n",
    "\n",
    "\n",
    "\n",
    "ax[1].fill_between(ej_th[:,0],commb_25,commb_35,color = 'grey',alpha = 0.4,\n",
    "             label = 'Best Fit: \\n 25 ppm - 35 ppm')\n",
    "\n",
    "\n",
    "\n",
    "#ax[1].plot(ej_th[:,0],commb_13,'k-',linewidth= 1,label = 'Minimum: 13 ppm'.format(ej_inc))\n",
    "\n",
    "\n",
    "ax[1].set_xlabel('Distance from Center of Imbrium', fontsize = 14)\n",
    "ax[1].set_xlim([500,5100])\n",
    "ax[1].legend(loc='upper right', fontsize = 14)\n",
    "\n",
    "\n",
    "ax[1].set_ylabel('Thorium (ppm)',fontsize = 14)\n",
    "ax[1].tick_params(axis='both', which='major', labelsize=12)\n",
    "\n",
    "\n",
    "\n",
    "ax[1].set_ylim([0,22])\n",
    "#ax[1].legend(loc='upper right',fontsize = 14)\n",
    "#ax[1].set_title('Minimum Concentration of Layer',fontsize = 14)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b939ae2f-e586-4a3f-b6c4-28efed781df2",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.savefig('imbrium_figure_25-35_2024.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9219aff0-6f02-4c21-93f8-bc87fbdf9104",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e74f929-0e21-4316-a642-7c4879708a98",
   "metadata": {},
   "outputs": [],
   "source": [
    "ej1 = kreep_val[np.argmin(sse_qs[0,:])]\n",
    "ej2 = kreep_val[np.argmin(sse_qs[1,:])]\n",
    "ej3 = kreep_val[np.argmin(sse_qs[2,:])]\n",
    "ej4 = kreep_val[np.argmin(sse_qs[3,:])]\n",
    "\n",
    "commb_no1 = (ej_th[:,2]*ej1 + ej_th[:,3]*0.63)/ej_th[:,1] \n",
    "commb_no2 = (ej_th[:,2]*ej2 + ej_th[:,3]*0.63)/ej_th[:,1] \n",
    "commb_no3 = (ej_th[:,2]*ej3 + ej_th[:,3]*0.63)/ej_th[:,1] \n",
    "commb_no4 = (ej_th[:,2]*ej4 + ej_th[:,3]*0.63)/ej_th[:,1] \n",
    "\n",
    "plt.plot(ej_th[:,0],commb_no1,'-',color = 'maroon', linewidth= 2,label = 'Q1 Best Fit: {:.1f} ppm'.format(ej1))\n",
    "plt.plot(ej_th[:,0],commb_no2,'-',color = 'orangered', linewidth= 2,label = 'Q2 Best Fit: {:.1f} ppm'.format(ej2))\n",
    "\n",
    "plt.plot(ej_th[:,0],commb_no3,'-',color = 'darkorange', linewidth= 2,label = 'Q3 Best Fit: {:.1f} ppm'.format(ej3))\n",
    "\n",
    "plt.plot(ej_th[:,0],commb_no4,'-',color = 'goldenrod', linewidth= 2,label = 'Q4 Best Fit: {:.1f} ppm'.format(ej4))\n",
    "plt.legend(loc = 'upper right')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42d10494-5805-449a-9dde-a4073ef6e91f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "607df3b9-c145-4a92-83e6-c29c4166d8ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "kreep_val = np.arange(1.,100.,0.5)\n",
    "\n",
    "front_bound = np.arange(0,3100,100)\n",
    "ejtotb = np.zeros_like(front_bound)\n",
    "ej_errsb = np.zeros((4,len(front_bound)))\n",
    "b = 0\n",
    "for bou in front_bound:\n",
    "    bounds = [bou,3300]\n",
    "    an_th = 0.63  ###from metzger 1973\n",
    "\n",
    "    sse_tot = np.zeros((1,len(kreep_val)))\n",
    "    sse_qs = np.zeros((4,len(kreep_val)))\n",
    "\n",
    "    idcs = np.where((dist_list < bounds[1])& (dist_list > bounds[0]))\n",
    "\n",
    "    for i in range(len(kreep_val)):\n",
    "        kr_th = kreep_val[i]\n",
    "\n",
    "        con = (ej_th[:,2]*kr_th + ej_th[:,3]*an_th)/ej_th[:,1] \n",
    "        s = interpolate.InterpolatedUnivariateSpline(ej_th[:,0], con) ###interpolate the data points with a univariate spline\n",
    "\n",
    "\n",
    "        er_qs = np.zeros((4,len(idcs))) + s(dist_list[idcs])\n",
    "\n",
    "        for q in range(4):\n",
    "            er_qs[q,:] += -thors[q,idcs][0] \n",
    "\n",
    "        sse_qs[:,i] = np.nansum(er_qs**2,axis = 1)\n",
    "\n",
    "        sse_tot[0,i] = np.nansum(er_qs**2)\n",
    "    min_err = np.min(sse_tot[0,:])\n",
    "\n",
    "\n",
    "    ejtotb[b] = kreep_val[np.argmin(sse_tot)]\n",
    "    ej_errsb[:2,b] = np.sort(kreep_val[np.argsort(np.abs(sse_tot[0,:] - 2*min_err))[:2]])\n",
    "    ej_errsb[2:,b] = np.sort(kreep_val[np.argsort(np.abs(sse_tot[0,:] - 4*min_err))[:2]])\n",
    "    b += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "835c3f2c-b5b6-4b53-ad95-d4fa066ad35c",
   "metadata": {},
   "outputs": [],
   "source": [
    "kreep_val = np.arange(1.,100.,0.5)\n",
    "\n",
    "back_bound = np.arange(1000,3100,100)\n",
    "ejtotf = np.zeros_like(back_bound)\n",
    "ej_errsf = np.zeros((4,len(back_bound)))\n",
    "b = 0\n",
    "for bou in back_bound:\n",
    "    bounds = [700,bou]\n",
    "    an_th = 0.63  ###from metzger 1973\n",
    "\n",
    "    sse_tot = np.zeros((1,len(kreep_val)))\n",
    "    sse_qs = np.zeros((4,len(kreep_val)))\n",
    "\n",
    "    idcs = np.where((dist_list < bounds[1])& (dist_list > bounds[0]))\n",
    "\n",
    "    for i in range(len(kreep_val)):\n",
    "        kr_th = kreep_val[i]\n",
    "\n",
    "        con = (ej_th[:,2]*kr_th + ej_th[:,3]*an_th)/ej_th[:,1] \n",
    "        s = interpolate.InterpolatedUnivariateSpline(ej_th[:,0], con) ###interpolate the data points with a univariate spline\n",
    "\n",
    "\n",
    "        er_qs = np.zeros((4,len(idcs))) + s(dist_list[idcs])\n",
    "\n",
    "        for q in range(4):\n",
    "            er_qs[q,:] += -thors[q,idcs][0] \n",
    "\n",
    "        sse_qs[:,i] = np.nansum(er_qs**2,axis = 1)\n",
    "\n",
    "        sse_tot[0,i] = np.nansum(er_qs**2)\n",
    "    min_err = np.min(sse_tot[0,:])\n",
    "\n",
    "\n",
    "    ejtotf[b] = kreep_val[np.argmin(sse_tot)]\n",
    "    ej_errsf[:2,b] = np.sort(kreep_val[np.argsort(np.abs(sse_tot[0,:] - 2*min_err))[:2]])\n",
    "    ej_errsf[2:,b] = np.sort(kreep_val[np.argsort(np.abs(sse_tot[0,:] - 4*min_err))[:2]])\n",
    "    b += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02ba5261-f8b8-43d8-89fb-6543128fc20c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize = [8,4])\n",
    "ax1 = fig.add_subplot(121)\n",
    "ax1.fill_between(front_bound,ej_errsb[2,:],ej_errsb[3,:],color = 'lightgrey',label = 'Error 200%')\n",
    "ax1.fill_between(front_bound,ej_errsb[0,:],ej_errsb[1,:],color = 'dimgrey',label = 'Error 100%')\n",
    "ax1.plot(front_bound,ejtotb,'k*',label = 'Best Fit')\n",
    "\n",
    "\n",
    "ax1.set_xlabel('Lower bound of best-fit region (km)')\n",
    "ax1.set_ylabel('$[Th_K]$ (ppm)')\n",
    "#plt.vlines(front_bound,ej_errs[2,:],ej_errs[3,:],color = 'green')\n",
    "ax1.set_ylim([10,55])\n",
    "\n",
    "ax2 = fig.add_subplot(122)\n",
    "\n",
    "\n",
    "ax2.fill_between(back_bound,ej_errsf[2,:],ej_errsf[3,:],color = 'lightgrey',label = 'Error 200%')\n",
    "ax2.fill_between(back_bound,ej_errsf[0,:],ej_errsf[1,:],color = 'dimgrey',label = 'Error 100%')\n",
    "ax2.plot(back_bound,ejtotf,'k*',label = 'Best Fit')\n",
    "\n",
    "\n",
    "ax2.set_xlabel('Upper bound of best-fit region (km)')\n",
    "#plt.vlines(front_bound,ej_errs[2,:],ej_errs[3,:],color = 'green')\n",
    "\n",
    "ax2.legend(loc = 'upper left')\n",
    "ax2.set_ylim([10,55])\n",
    "ax2.set_yticklabels([])\n",
    "fig.subplots_adjust(wspace = 0.05)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86e7e57f-82dc-452b-8769-ca7798b0f84b",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.savefig('/Users/janie/Documents/presentations/varying_best_fit_bounds.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12f3f4c9-7b68-4c60-aac2-493a78e2c2e5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d51198c7-659d-4d7b-844c-5f521824b70d",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sort([4,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbdbcb47-cd8c-40de-b659-aa0078880adf",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "plt.fill_between(back_bound,ej_errs[2,:],ej_errs[3,:],color = 'lightgrey',label = 'Error 200%')\n",
    "plt.fill_between(back_bound,ej_errs[0,:],ej_errs[1,:],color = 'dimgrey',label = 'Error 100%')\n",
    "plt.plot(back_bound,ejtot,'k*',label = 'Best Fit')\n",
    "\n",
    "\n",
    "plt.xlabel('Upper bound of best-fit region (km)')\n",
    "plt.ylabel('$[Th_K]$ (ppm)')\n",
    "#plt.vlines(front_bound,ej_errs[2,:],ej_errs[3,:],color = 'green')\n",
    "\n",
    "plt.legend(loc = 'upper left')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ac80af6-6343-4f0a-afd7-150125449b93",
   "metadata": {},
   "outputs": [],
   "source": [
    "kreep_val = np.arange(1.,100.,0.5)\n",
    "\n",
    "front_bound = np.arange(0,3100,100)\n",
    "back_bound = np.arange(2000,3300,100)\n",
    "\n",
    "fbb,bbb = np.meshgrid(front_bound,back_bound)\n",
    "\n",
    "ejtot = np.zeros((len(front_bound),len(back_bound)))\n",
    "b = 0\n",
    "for f in range(len(front_bound)):\n",
    "    for b in range(len(back_bound)):\n",
    "        if back_bound[b]-front_bound[f] > 500:\n",
    "            bounds = [front_bound[f],back_bound[b]]\n",
    "            an_th = 0.63  ###from metzger 1973\n",
    "\n",
    "            sse_tot = np.zeros((1,len(kreep_val)))\n",
    "            sse_qs = np.zeros((4,len(kreep_val)))\n",
    "\n",
    "            idcs = np.where((dist_list < bounds[1])& (dist_list > bounds[0]))\n",
    "\n",
    "            for i in range(len(kreep_val)):\n",
    "                kr_th = kreep_val[i]\n",
    "\n",
    "                con = (ej_th[:,2]*kr_th + ej_th[:,3]*an_th)/ej_th[:,1] \n",
    "                s = interpolate.InterpolatedUnivariateSpline(ej_th[:,0], con) ###interpolate the data points with a univariate spline\n",
    "\n",
    "\n",
    "                er_qs = np.zeros((4,len(idcs))) + s(dist_list[idcs])\n",
    "\n",
    "                for q in range(4):\n",
    "                    er_qs[q,:] += -thors[q,idcs][0] \n",
    "\n",
    "                sse_qs[:,i] = np.nansum(er_qs**2,axis = 1)\n",
    "\n",
    "                sse_tot[0,i] = np.nansum(er_qs**2)\n",
    "            min_err = np.min(sse_tot[0,:])\n",
    "\n",
    "            ejtot[f,b] = kreep_val[np.argmin(sse_tot)]\n",
    "        else:\n",
    "            ejtot[f,b] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "864feeb2-de5e-46eb-99be-2627544fcc05",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = plt.pcolormesh(back_bound,front_bound,ejtot)\n",
    "\n",
    "plt.colorbar(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "348c62d9-d620-4020-b16c-e1af09d798e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "####lower bound on concentration of imbrium reservoir \n",
    "lb_th = 13. ###lower bound on thorium\n",
    "### concentration of sub-imbrium reservoir \n",
    "an_th = 0.63\n",
    "r0_imb = 360.5 ###transient radius\n",
    "h_exc_imb = r0_imb/5\n",
    "\n",
    "tot_vi = 0.1*np.pi*r0_imb**3 ### simplification - not worrying bout the 10% falling back in: assume well mixed\n",
    "\n",
    "x_imb = 28. #thickness of crust (Miljkovic 2016)\n",
    "x_pm = 7. ##plus or minus\n",
    "\n",
    "hk_imb = h_exc_imb - x_imb ##pm x_pm (7)  \n",
    "\n",
    "\n",
    "\n",
    "###radius of bottom paraboloid\n",
    "rk_imb = r0_imb * np.sqrt(hk_imb)/(np.sqrt(h_exc_imb)) ### pm x_pm/2\n",
    "\n",
    "####error propagation\n",
    "\n",
    "hk_rse = x_pm/hk_imb ##relative error\n",
    "\n",
    "rk_pm = r0_imb*(hk_imb* hk_rse/2)/np.sqrt(h_exc_imb) ## square root of number means divide pm by relative error\n",
    "\n",
    "rk_pm2 =  rk_imb**2*((rk_pm/rk_imb)/2) ### error of rk_imb**2\n",
    "\n",
    "\n",
    "##### volume of bottom paraboloid\n",
    "\n",
    "v_kreep_imb = np.pi*(rk_imb**2)*hk_imb/2\n",
    "\n",
    "v_k_pm = np.pi*rk_imb**2*hk_imb*np.sqrt((rk_pm2/rk_imb**2)**2 + hk_rse**2)/2\n",
    "\n",
    "#### volume of crust above\n",
    "v_crust = tot_vi - v_kreep_imb \n",
    "\n",
    "\n",
    "#### concentration of bottom reservoir\n",
    "th_k_lb = (tot_vi*lb_th - v_crust*an_th)/v_kreep_imb \n",
    "\n",
    "#th_imb = (tot_vi*no_ej + v_crust * an_th)/v_kreep_imb *np.sqrt(((v_k_pm*an_th)/(tot_vi*no_ej + v_crust * an_th))**2 + (v_k_pm/v_kreep_imb)**2)\n",
    "\n",
    "th_k_pm = (tot_vi*lb_th -v_crust * an_th)/v_kreep_imb *np.sqrt(((v_k_pm*an_th)/(tot_vi*lb_th - v_crust * an_th))**2 + (v_k_pm/v_kreep_imb)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83bf9137-e67f-4e2d-b39f-711d9519f14a",
   "metadata": {},
   "outputs": [],
   "source": [
    "th_k_lb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "025461ae-cad8-4360-99fe-c6542c31e45a",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(v_kreep_imb, '$\\pm$', v_k_pm)\n",
    "\n",
    "print(th_k_lb, r'$\\pm$', th_k_pm) ####parts per million\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
